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1. Introduction 

Simulations of charged systems face a big computational challenge due to the long- 
range nature of the electrostatic interaction. If N is the number of charges, then the 
computational cost of the most naive approach to evaluate the interaction energy would 
scale as N"^, since every charge interacts with every other charge. Very sophisticated 
algorithms have been developed to tackle this problem and to reduce the computational 
complexity. The most prominent ones are the so-called P^M method ( "particle-particle 
/ particle-mesh" ) , which is based on Fast Fourier Transforms and scales as A^ log A^ PP , 
and the Fast Multipole method |2I which scales linearly with A^. 

A similar problem arises in the simulation of Brownian particles which interact 
hydrodynamically: Their stochastic displacements are highly correlated, due to fast 
diffusive momentum transport through the solvent. For sufficiently slow particles, a 
quasi-static approximation works excellently, and in this case the correlation function 
decays as 1/r (r interparticle distance) |3j, just as in electrostatics. For these systems, it 
has turned out that it is both much simpler and also more efficient to explicitly simulate 
the momentum transfer through the surrounding solvent. This makes the simulation 
of several ten thousands of Brownian particles feasible jH Ej- Although most of the 
computational effort goes into the flow field (for two reasons — one needs reasonable 
spatial resolution of the flow fleld, and it moves much faster than the Brownian particles), 
this approach ultimately wins, because it is inherently local, and therefore scales linearly 
with A^. 

This observation raises the question if something similar could be tried for Coulomb 
interactions. After all, electrostatics is just the quasi-static limit of full electrodynamics. 
The obvious approach would be to couple a system of charges to an electromagnetic fleld 
which propagates according to the Maxwell equations (ME), and then run Molecular 
Dynamics (MD). A suitable acronym for such a method might be MEMD ("Maxwell 
equations Molecular Dynamics"). Just as in the hydrodynamic case, this is an 
intrinsically local algorithm, and therefore scales linearly. The instantaneous 1/r 
interaction is thus replaced by some retarded interaction travelling with the speed of 
light c. Using the actual physical value of c will of course not work, since then the 
separation of time scales between charges and flelds will be prohibitive. However, there 
is no need to take such a large c value. It is sufficient to just make c large enough 
such that the quasi-static approximation still holds to sufficient accuracy. This is the 
lesson we have learned from Car-Parrinello (CP) simulations jSj, where the electrons are 
assigned an unphysically large mass, precisely for the same reason. The analogy between 
MEMD and CP actually goes much further, as we will see below. This should not be 
too much of a surprise, since the universal applicability of the CP approach to a wide 
variety of problems in physics (e. g. classical field theories) has already been observed in 
the original publication jB] , and exploited in the context of classical density-functional 
theory [Zj. 

The MEMD idea has been pursued recently by A. C. Maggs and collaborators 
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[HI El UD] , and by us, in close contact with him. He has made a couple of very important 
observations, which have deepened our insight into the approach significantly, and 
contributed to the answer of a number of very important questions: 

(i) Is Maxwell dynamics the only possible way to propagate the fields? The answer 
is no; it is also possible to propagate them in a diffusive fashion. This has been 
implemented by means of a Monte Carlo algorithm [Hlini for a lattice gas of charges. 

(ii) If we restrict attention to Hamiltonian or quasi-Hamiltonian dynamics of the 
system, and want wave-like propagation of the signal, is then Maxwell-style 
dynamics the only choice? The answer is a cautious yes; one can show that the 
Maxwell equations arise in a very natural way if one derives the method along the 
lines of CP. 

(iii) Is there a contradiction between the Lorentz covariance of the ME, and the strictly 
nonrelativistic setup of MD? The answer is no; the Lorentz covariance actually has 
to do with the fact that the value of c is the same in all reference frames. This, 
however, is not the case here: In our context, c means nothing but the propagation 
velocity of electromagnetic waves relative to the discretization lattice which provides 
an absolute reference frame (an "ether"). 

(iv) Is it necessary to use a large value of c to avoid violation of a quasi-static 
behavior? The answer is no as long as just static properties of the system in 
thermal equilibrium are considered — the values of these properties turn out to be 
completely independent of c. 

(v) Is it necessary to apply a thermostat to the system? Ref. (TUI claims yes, in order 
to avoid unwanted conserved quantities. Our belief is no, based upon the fact that 
the particle dynamics provides lots of nonlinearities into the equations of motion. 
For more details, see below. 

(vi) How is MEMD implemented? The previous papers have been rather brief on this 
issue; we try to provide somewhat more detail. 

(vii) How does MEMD perform, in particular in comparison with existing methods? In 
this respect, there is also so far only little information available. In the present 
paper, we report some benchmark results which give us some feeling for the quality 
of the approach — although these are quite preliminary, and still far from providing 
a clear and comprehensive picture. 

In what follows, we will essentially re-derive the MEMD algorithm put forward in Ref. 
|1UJ . and discuss some details of our implementation, which differs slightly from that of 
Ref. ^01 • We will then present some benchmark results, comparing MEMD with P^M for 
the same system. For our chosen set of parameters, we find rather similar or even better 
computational efficiency. However, this comparison should not yet be considered as the 
final answer: Firstly, MEMD can probably still be speeded up significantly by combining 
it with a direct evaluation of Yukawa-like forces on short length scales, roughly along 
the lines as suggested in Ref. [10]. Secondly, the dependence on the thermodynamic 
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state point (in particular, on density) has not yet been investigated. Physically, it is 
clear that the efficiency of MEMD depends (i) on the number of operations required 
to propagate the system for one time step, and (ii) on the time needed to build up 
electrostatic correlations on the relevant length scale, which is, in essence, the Debye 
screening length. For this reason, one should expect that the efficiency depends rather 
strongly on the speed of light, and also on the density, since the Debye length decreases 
as a function of density. In other words, one expects MEMD to work particularly well 
for rather dense systems simulated with a large c value — and this is precisely the 
regime where our preliminary comparison has been done. 

2. Continuum Theory 

We start out from Maxwell's equations in vacuum, using standard SI units: 

V-E =-p (1) 

V-H =0 (3) 

-^ ^ ^ dE 

VxH = j + e^ — , (4) 

— * — * 

where eo is the vacuum dielectric constant, c the speed of light, E the electric field, H 
the magnetic field, p the charge density and j the current density, which are coupled 
via the continuity equation 

l + V.J^O, (5) 

Electrostatics is obtained by setting the current and all time derivatives to zero, implying 
that the magnetic field vanishes: 

V -E =-p (6) 

V X E = 0. (7) 

Note that this set of equations also results from taking the limit c ^ oo. This means 
that the electric field will be just an electrostatic field as long as the charges move at 
much slower velocity than c. Furthermore, the Lorentz force on a charge e, 

FL = e(E + -^v X H^ (8) 

[v denoting the charge's velocity) will just reduce to the electrostatic force eE in the 
same limit. This is the justification of the fact that MD simulations of common materials 
usually treat interactions of charges just as electrostatics. In turn this means that one 
will obtain electrostatic behavior whenever c is large compared to all particle velocities, 
as already stated in Sec[T] 
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The conventional approach tries to find a solution to Eqs. iniand[7| for a given charge 
distribution p. Note that both equations must be satisfied for strict electrostatics, since 
Eq. ini only fixes the longitudinal component of the electric field, while the condition of 
vanishing transversal component is coded in Eq. [7| No local way of finding the solution 
directly is known. 

As a first step, we re-formulate the electrostatic problem in terms of a constrained 
variational problem. Gauss' law (Eq. ^ is viewed as a constraint which selects a certain 
surface out of the space of electric field configurations; we will call this the "constraint 
surface" (CS). We now minimize the electric field energy, 

^0 f .3 - B2 



'Hef = -^J d'rE' (9) 

under the constraint Eq. El This can be done as follows: Suppose Eo is some field on 
the CS, where a non- vanishing transversal component is admitted. Then all fields on 
the CS can be written in the form 

E = Eo + V xQ, (10) 

— * 

where is allowed to pass through all field configurations without any restriction. We 
thus write T-Cef in terms of the B field, 

Hef ({e}) = I / d'r(Eo + V X e)' (11) 

and the minimum condition as 

^n,, ({§}) = (12) 

or 

= V X (^0 + V X e) = V X E, (13) 

i. e. Eq. |7| The variational problem is thus seen to be equivalent to the original 
electrostatic problem. We can say that the system is on its Born-Oppenheimer surface 
(BOS) if, for given charge distribution p, the electric fields are on the constraint surface, 
and the field energy is minimal. 

Alternatively, one may also look at the problem in Fourier space: Let the Fourier 
transform of E{f) be defined as 

E{k) = (27r)~^/2 f d^rEif) exp (-ik ■ f) (14) 

and let k denote the unit vector in the direction of k. Then we can decompose the electric 
field into a longitudinal component E\\ and a transversal component E^, E = E\\+ E± 
with E\\ ■ k = E\\ and E± ■ k = 0. Then Eqs. iniand[7|are transformed to ik ■ E = p/tQ 

-- 0. Furthermore, the electric field energy 



and ik X E = Q, OT E\\ = p/{ik 
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Again, one sees that the longitudinal component is determined by the charge 
distribution, while the transversal component is just minimized away. 

In the spirit of CP simulations, we now wish to replace the precise solution of the 
minimization problem by some (to a certain degree artificial) dynamics which keeps 
the system precisely on the CS but allows fluctuations around the BOS. An interesting 
observation by Maggs [H] is that arbitrarily large deviations from the BOS are permitted 
as long as one does statistical mechanics in the canonical ensemble, and is interested 
in static properties only. This is easily understood by looking at Eq. ^1 There one 
sees that the total Hamiltonian decomposes into two additive contributions: The first 
term 7i\\ is just the standard electrostatic Coulomb Hamiltonian, while the second term 
7i± is the energy stored in the additional transversal degree of freedom, describing the 
amount of deviation from the BOS. Additivity, however, implies that the Boltzmann 
factor factorizes, 

exp i-m) = exp {-mw) exp i-l3n±) (16) 

[jS = l/{kBT), where T is the absolute temperature and fc^ Boltzmann's constant), 
which in turn means that the (artificial) transversal degree of freedom is statistically 
independent from the physical longitudinal one, and hence does not affect statistical 
averages of observables which only depend on the charge configuration. This will be 
worked out in some more detail at the end of this section. 

Having relaxed the condition Eq. [3 we now turn our attention to Eq. IHl Suppose 
that at time t = we have found the full solution to Eqs. IHl and [7| by some (slow) 
procedure; we call this solution Eo(t = 0). This is obviously on the CS. The system will 
then stay on the CS if the time derivative of Gauss' law vanishes: 

V-^--p = 0. (17) 

Now as the dynamics proceeds, the continuity equation, Eq. El will automatically hold 
as long as charges are moved around in the simulation cell by a local updating scheme. 
This allows us to re-write Eq. ITTIas 

V-(^ + ^j)=0. (18) 

We can therefore use the current density to straightforwardly construct an electric field 
which stays on the CS. One just has to integrate E = —j/tQ in time; this is a manifestly 
local updating scheme. We thus obtain 

Eo(t) = E^{t = 0) - - fdr]{r). (19) 

This solution obviously is on the CS, but unfortunately not the correct electric field 
(example: For a constant ring current with vanishing charge density one would obtain 
an electric field which grows linearly in time). We therefore generalize this, as before, 

to 

i(t) = io(t) + Vx e(t) (20) 
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with B(t = 0) = 0. However, we now do not minimize Ti-EF with respect to G, but 
rather supply some artificial dynamics to this field. There is no unique way of doing 
this. One possibility is to postulate an overdamped relaxational dynamics governed 
by Ti.EF] this has been explored in detail in Ref. [Hj. In the present paper, we rather 
study, as in Ref. [10], a CP-style dynamics, where the equation of motion for 9 is of 
second order in time. We thus need to supply an initial condition for G, too; we choose 
0(t = 0) = 0. The most straightforward way to generate a coupled dynamics is to add 

a kinetic energy term (l/2)(eo/c^) J d^fQ to the system Lagrangian; here the prefactor 
is a mass-like parameter, to be freely chosen in analogy to the electron mass in CR c 
will later on turn out to be the speed of light. 

Since Eq depends on the charge distribution in a not very straightforward way, it 
is more convenient to rather write the Lagrangian in terms of the total field E, and to 
take into account the integration of Eq = — j/eo by means of a non-holonomic constraint 
which keeps the system on the CS: 

eoE = eoEo + eoV x 9 = -J + cqV x 9, (21) 



1. e. 



eo^ + /-eoVx 9 = 0. (22) 

This is nothing but the fourth Maxwell equation, Eq. 01 if we identify H = eo9. We 
thus see that the continuity equation (Eq. EJ as well as the first and fourth Maxwell 
equation (Eqs. [T] and E} are built into the scheme regardless of the details of the 9 
dynamics. 

Denoting the particle masses with nii, their coordinates with r^, and the 
interparticle potential (of non-electromagnetic type) with U, we can thus write the 
Lagrangian as 

L=T.Y''i'-U (23) 



+ II. I d^fQ - 



2c2 y 2 7 

+ Jd^rA-(^eoE + j-eoVxQy, 

here the field A is a Lagrange multiplier; it will later on turn out to be the vector 
potential. Such a constrained variational problem with Lagrange multipliers can 
be treated by Arnold's so-called "vakonomic" ("variational of the axiomatic kind") 
formalism JT]. The recipe how to do variational calculus within that formalism is very 
simple: One just has to treat all occuring variables, including the Lagrange multipliers, 
as if they were independent degrees of freedom. It is thus straightforward to obtain 
the equations of motion. Variation with respect to A just yields the fourth Maxwell 
equation, Eq. |221 Variation with respect to E yields 

E = -A, (24) 
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while from variation with respect to G we obtain 

\e = V xA. (25) 

This is equivalent to the remaining two Maxwell equations, Eqs. El and 01 Inserting Eq. 
EH plus H = eoO, into Eq. (23 we obtain directly Eq. |21 Furthermore, we can integrate 
Eq. |2niin time, which, together with H = eoQ, and the initial condition {A and G both 
vanish at time t = 0) yields 

H = eoc^V X A. (26) 

Taking the divergence of this equation, one directly obtains Eq. JHl The interpretation 
of A as the vector potential is also obvious from Eqs. j^ and |2S1 since these are the 
standard relations between the electromagnetic fields and the vector potential. It should 
be noted that our derivation has led us in a natural way to the so-called temporal gauge 
P^ where the scalar potential vanishes identically, and there is no restriction on A. 

For deriving the equations of motion for the particles, we first note that charge and 
current densities are written as 

p=Y,e,6{f-f,) (27) 

i 

j = ^^eirid^f-ri) (28) 

i 

where Cj is the charge of the ith particle. Hence the current term in the Lagrangian is 
written as 

l£fA-j=J2e,A{n)-r',. (29) 

i 

After a few lines of algebra one then finds the particle equations of motion: 

mih = -% + h, (30) 

where the Lorentz force F^ is given by Eq. jHl 

To summarize: The requirement of local updates, combined with treating the 
deviations from the BOS in the CP manner, has led us in a natural way to standard 
electromagnetism, where the temporal gauge turns out to be the most appropriate one 
for our purposes. It should be stressed that this is a consistent non-relativistic setting, 
where the equations of motion are valid in one particular chosen frame of reference. 

As it is common practice in electromagnetism ^3], we can now simplify the 
Lagrangian treatment by considering A as the (only) field degree of freedom, while 
E and H are derived quantities according to Eqs. Ol and EEl The dynamical system 
of charges and electromagnetic field is then completely described by (i) the equation of 
motion for the particles, Eq. (201 and (ii) the fourth Maxwell equation, Eq. |3J which is 
the inhomogeneous wave equation for A: 

-A = -c^V X (V X /) + -J. (31) 



dt^ V y ' eo 
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To derive these two equations of motion, it is sufficient to consider the Lagrangian 

^ = E f -"' - ^ (32) 

i 



2 



fd^f(yx ly + - fd^fA 



+ jd^rA-j. 



This dynamics has a couple of very desirable properties: Firstly, since the dynamics is 
manifestly Hamiltonian (it is derived from a Lagrangian), it conserves the phase-space 
volume and the energy, the latter being given by 

'>^=T.Y^^+U (33) 



2 J 2eoc2 J 



Furthermore, one can show that the total momentum, given by 

P = J2 ^ih + — f d^rE X H, (34) 

is conserved as well. For the proof one can employ the dynamic equations for the 
particles and fields, and make use of the identity 

f d^fX X (V xX) = f d^fx{V-X), (35) 

which holds for any vector field X as long as partial integration with vanishing boundary 
terms can be applied. 

At this point, we modify the equations of motion by discarding the magnetic force 
on the particles, 

dU 

miTi = -iT^ + eiE{r.i). (36) 

This simplifies the algorithm significantly, while the most important features still hold. 
Of course, this modified dynamics is no longer Hamiltonian. Nevertheless, the energy, 
as given by Eq. EHl is still conserved. Furthermore, the (properly defined) phase space 
volume is also conserved. In order to see this, we first write the equations of motion in 
pseudo-Hamiltonian style as 

— r, = —Pi (37) 

dt rui 

|p. = - II + e,i(rl) (38) 

p^-E (39) 

^E = c^V X (V X 1) - Ij, (40) 

where the Pi are the kinematic (and not the canonically conjugate!) particle momenta, 
and the fields A and E (roughly) play the roles of coordinates and momenta, respectively. 
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Now, phase space volume conservation for some (high-dimensional) dynamical 
system, given by the equation 

f=/(f), (41) 

where x comprises the set of all phase-space variables, holds if and only if V • / = 
(in analogy to incompressible flow in hydrodynamics). It is trivially checked that this 
relation does hold for our system. 

However, momentum conservation does not hold for our modified dynamics. The 
momentum carried away by the electromagnetic waves is not completely balanced by 
the particle momenta. Rather, we have the relation 

Y.Pi = const. + O (c"2^ . (42) 

i 

This is not a catastrophe, since momentum conservation is usually only important in 
studies of dynamics. However, for such calculations one has to use a fairly large value 
of c anyways, since otherwise the electromagnetic field is not in its quasi-static limit, 
and the particle trajectories get too much distorted. Furthermore, one must expect that 
momentum conservation is also violated as a result of the lattice discretization, which 
breaks the translational invariance of the system. 

We now assume that the dynamics is sufficiently nonlinear to make the system 
ergodic. This seems reasonable for the case of a many-charge system, in particular if 
the potential U has a strongly repulsive core to facilitate "collisions". We therefore 
assume that the system has no further important conservation law except for the fact 
that it stays on the CS, and that the energy H is conserved. The additional conserved 
quantities mentioned in Ref. JU] probably apply only to the charge-free case, in which 
the system is harmonic and hence integrable. We can hence apply standard statistical 
physics to the system and assume that the dynamics results in an equidistribution of 
states in terms of the variables fi,pi,A and E (microcanonical ensemble). Making use 
of the fact that thermodynamic ensembles are equivalent in the large-system limit, we 
can instead employ the canonical ensemble, which is easier. With (3 = l/{kBT), where 
ks is Boltzmann's constant and T the absolute temperature, we may therefore write 
the partition function as 



dfi f dpi jvAJVE exp {-(371) 



x<5(v-E--t-p), (43) 

where Ti is given by Eq. EHl It is now straightforward to integrate out the momenta, 
the A field, and the transversal component of the E field. The integration over 
the longitudinal component of E cancels with the delta function, such that the only 
remaining degrees of freedom are the particle coordinates, for whose potential of mean 
force we hence find 

n,onf = U+'-^Jd^fE^; (44) 



n^onf = U + -,j^ld^r I d\- l^i^. (45) 
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here E is nothing but the solution of the standard electrostatic problem, Eqs. IHl and [7| 
i. e. the Coulomb field. Inserting this field into Eq. EJ we find the standard Coulomb 
Hamiltonian, 

2ATTeoJ J If-P 

This demonstrates that the particles behave statistically in the same way as if they would 
directly interact Coulombically. This concludes the derivation. On the lattice, however, 
we have to take into account that the above Hamiltonian includes unphysical self- 
interactions, which we have to subtract (without such a subtraction the self-interaction 
would ultimately, i. e. in the continuum limit of vanishing lattice spacing, completely 
dominate the behavior), and that instead of the 1/r^ Coulomb field we have to insert 
the lattice-discretized solution of the lattice equations corresponding to Eqs. El and [3 
This shall be discussed in the next section. 

3. Discretization, Lattice Green's Function and Self— Interaction 

For implementation on the computer, the equations need to be discretized with respect 
to both space and time. For the moment, we will only consider the spatial discretization, 
and consider time still as a continuous variable. The issue of time discretization is 
discussed in [Appendix B 



We use a spatial discretization scheme jHl IHI where the charges are interpolated 
linearly to the eight surrounding lattice sites of a simple-cubic lattice. The currents, as 
well as the fields E and A are put onto the connecting links. The curl of link variables 
is put onto the lattice plaquettes, and the curl of plaquette variables onto the links (in 
both cases one uses the four fields which encircle the result). Furthermore the divergence 
of link variables is put onto the sites, using the adjacent fields, while the gradient of a 
scalar variable located on the sites is a link variable. For more details, see [Appendix A| 

Let us now discuss how the Coulomb potential looks on the lattice. Obviously, we 
have to solve Eqs. [Bland [3 on the lattice. As in the continuum, we can take into account 
the longitudinal character of the electric field by the ansatz 

E = -V0, (46) 

where is the electrostatic potential on the sites, and V is the lattice-discretized 
gradient. We thus obtain the Poisson equation on the lattice, 

- VV = -P, (47) 

Co 

where the lattice version of the operator V^ is clear from the previous definitions of 
gradient and divergence. 

A system with periodic boundary conditions is invariant with respect to lattice 
translations, and this allows us to write 

0(f) = — ^ G(f - f')p{f^), (48) 

^0 ^ 
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where a is the lattice spacing, r denotes the sites, and G is the lattice Green's function, 
obeying the equation 

-V'G{r) = \6{r), (49) 

where (5(f^ is the Kronecker symbol, i. e. S{f) = 1 for r = 0, and 6{f) = for all other 
lattice sites. 

For an L^ x Ly x Lz lattice with periodic boundary conditions the solution 
can be obtained straightforwardly via discrete Fourier transformation. At the site 
r = a{nx,ny,nz) one finds 

Lx-l' 



a (r) = £ exp (2vrz^) 



Px=0 

Ly-l' / 

X ^2 6xp 27ri 






Py=o \ ^y / 

X }_^ exp 27r2-— 

X . . J G{px,py,pz), (50) 

where J2' indicates that {px,Py,Pz) = (0,0,0) is excluded (for reasons of overall charge 
neutrality), and G is given by 



G{px,Py,Pz) ^ = 6- 2cos f27r-pj 



-2cos(27ry^J (51) 

o fo P^ 

— 2 cos [ ZTT 

V Lz 

A lot is known about this function, in particular in the limit Lj — i> oo fHl QEl EI- 
For our purposes, however, it is sufficient to note that (i) G can be calculated at the 
beginning of the simulation once and for all, including the finite size effect, and that (ii) 
G{r = 0) is finite, even in the hmit Lj — > cxd (but, of course, keeping a fixed). 
We thus find for the potential of mean force (cf. Eq. 



nconf = U + l—Y.T.G{r- f')q{r)q{n, (52) 

where g(r) = a^p{r) is the charge on site r. Now, the charges on the sites are related 
to the charges Cj on the particles via the interpolation scheme, g(r) = I]iejs(r, Tj) and 
qir') = J2j^js{'r',rj), where s is the "smearing" function. Inserting this into Eq. |S21 
we find an effective interaction between different particles i y^ j, but also an unphysical 
self-energy term for i = j. This is given by 

Useif, = ^— EEG'(r -r-)e^(r,rl).(r-,rl). (53) 
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This depends explicitly on the particle coordinate Tj. The physical interpretation is 
simply that the Coulomb repulsion from the interpolated charges on the cube corners 
tries to drive the particle into the center of the cube. For small lattice spacings, this 
effect dominates over all other interactions, and therefore must be taken care of. For 
this reason, we add a term — J2i Useif,i to the interparticle potential U , and apply the 
corresponding force to the particles. This is feasible since G is known explicitly, and the 
smearing function s is short-ranged, such that X^r runs over eight sites only. 

4. Yukawa Subtraction 

Rottler and Maggs ^U] suggest another subtraction scheme which has the nice 
property of introducing another optimization parameter k into the method. Essentially, 
interactions up to the length scale k~^ are done in real space, while only the residual 
long-range part beyond k~^ is treated via the dynamics. The disadvantage, however, 
is that it does not treat the lattice effects completely rigorously. We hence believe that 
probably the best method consists of a combination between our lattice Green's function 
subtraction, and their "dynamic Yukawa" approach. 

In order to understand the latter, let us first consider the functional 

2 
and study, for fixed p, 

^-°- (-) 

It is straightforward to see that (i) this variational problem is equivalent to the Poisson 
equation for the electrostatic potential 0, and that (ii) insertion of the solution into 
JF yields JF = +(1/2) /(i^rp0, i. e. the correct electrostatic energy. However, this 
functional is useless for dynamic simulations where one would try to simulate a coupled 
dynamics of p and 0. The reason is that the V0 term has the wrong sign, such that 
arbitrarily large variations of are favored and the simulation would be inherently 
unstable (the partition function for integrating out the field would not exist). 

— * 

A well-behaved theory, however, is obtained by just turning the sign of the V0 
term: 

fO [ _,3_-> /^ i\2 

2 

This results in +V^0 = p/eo, and insertion into the functional yields again JF = 
+ (1/2) J d^fp(j). Since, however, is just the negative of the real (physical) electrostatic 
potential, one obtains a theory which describes attraction between like charges and 
repulsion between unlike charges. We now introduce an additional field degree of 
freedom 0, and couple this to the original method (Lagrangian) via 

^0 f ,3-i2 



J^=-- fd^f{V(j)y + fd^fp(j) (54) 



j^=+^A j d^r{y(fy + j d^^pcj), (56) 



L + i^ I d'fr (57) 



fo 
2 



I d^r {Vcj))^ - f d^f, 
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Here c^p is another dynamical parameter of dimension velocity. It can be set identical to 
c, but need not. This modified method would result in an additional potential of mean 
force between the charges which would exactly cancel the original Coulomb potential 
(including self-terms). This is apparently not useful. However, we can introduce a 
slightly modified coupling with a screening parameter k > 0: 

L^L + ^ f£r^^ (58) 

2c^ J 

^0 /"j3-?/^,\2 ^0 f ,3-f,2i2 



— / d'f {y^j -— d^fK^cj)^ - / d^rpcj). 

This introduces an additional potential of mean force between the charges, which, in 

the continuum limit, would read 

1 C'e ' 

Uyir) = ^^exp(— Kr), (59) 

^ ^ Aireo r ^ ' ^ ' 

such that unlike charges repel each other with a screened Coulomb interaction. This 

weakens the original Coulomb interactions on a local scale, and can be corrected by 

adding —Uy to the standard interparticle potential. Here one can use the continuum 

version of the potential; this will only serve to decrease the infiuence of lattice artifacts. 

In principle, this also alleviates the self-energy problem. However, the lattice 

Green's functions of the unscreened and screened Coulomb case are slightly different 

and only coincide in the limit k — ;► 0. In preliminary tests we found that fairly 

small screening parameters are needed to overcome the self-energy problem with high 

accuracy. We therefore believe that one should rather try to subtract the self-energy 

for both the unscreened and the screened interaction separately by the respective exact 

lattice Green's function. In this case, the Yukawa subtraction would no longer serve 

the purpose of overcoming self-energies, but rather to resolve interparticle interactions 

rather faithfully on a local scale, such that (hopefully) larger lattice spacings are feasible. 

Further investigations are necessary on this issue. 

5. Numerical Results 

As a simple test system, we have studied A^ charged particles in a cubic box with periodic 
boundary conditions. They interact via a purely repulsive Lennard- Jones (LJ) potential 



U 



LJ 



4e 


\r ) 


\rj 


6 1 

+ 4 


r < l^l^o 











r > 2^/^a 



(60) 



We choose a unit system where the potential parameters a and e, as well as the particle 
mass m, are set to unity. Time is thus measured in units of t^j = Jma'^/e. We 
study systems at temperature ksT = 1 and particle number density p = 0.07. The 
equations of motion were integrated by the algorithm outlined in [Appendix B| (no 
Yukawa subtraction), using a time step h = 0.01. The friction constant for the Langevin 
thermostat was set to C = 1. 
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Figure 1. Pair correlation function of like charges at density p = 0.07 and Bjerrum 
length Ib = 20, comparing data obtained with P'^M with those from MEMD for 
different lattice spacings. 



Each particle is assigned a charge ±e. The strength of the electrostatic interaction 
is given in terms of the Bjerrum length 

^« = r4r^- ^61) 

We first started out with the value Ib = 20 (rather strong electrostatic couphng). We 
chose this system because it had been studied previously by P^M ^H]- However, it has 
turned out that this is not the best state point for a benchmark, since the coupling is 
so strong that it actually induces phase separation (gas-liquid transition). This is in 
accord with the phase diagram presented in Ref. ^H]; the system studied there is not 
too different from ours. 

The number of particles was set to A^ = 2000. Both the P^M and the MEMD 
calculations were done within the framework of the "ESPResSo" software package ^H] 
of the Theory Group at the MPI for Polymer Research, Mainz. In both cases we used a 
program version which was fully parallelized, based upon domain decomposition. The 
P^M parameters were optimized using an automatized routine building upon the work of 
Refs. pUl l21j. where a formula for the relative error of the force per particle, AF/F, was 
derived. The routine provides optimized simulation parameters after an upper bound for 
AF/F has been supplied. For our system, we required an accuracy of 10~^, resulting in 
the following parameters: Mesh size 32^; 5th order charge assignment; real-space cutoff 
8.2; a = 0.36 (this parameter controls the split-up of the computational load between 
real and Fourier space). For the MEMD calculations, we used c = 1 and varied the 
lattice spacing a. 
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Figure 2. Pair correlation function of unlike charges at density p = 0.07 and Bjerrum 
length Ib = 20, again comparing P^M with MEMD. 



The pair correlation functions of this system are shown in Figs. Q and El The runs 
were long enough to equilibrate the system reasonably well on the local scale. As a 
control, we also ran a more accurate P^M simulation and found no visible difference 
from the original P^M result. The MEMD results in turn confirm that the static 
correlations converge towards those of the real electrostatic system when the lattice 
spacing a decreases. We found a value of a = 0.53 acceptable, corresponding to a 58^ 
lattice for the N = 2000 system. For such a fine lattice, there is practically never more 
than one particle per cube. 

We did not study this system further, since benchmarks at this state point are 
severely hampered by the gas-liquid transition: On the one hand, the system needs 
a long time to equilibrate (i. e. to condense macroscopically) , and on the other hand 
the particle density is very inhomogeneous in the relaxed state. This, in turn, is very 
detrimental to efficient geometric parallelization, since some processors have to treat 
very many particles, while others work on essentially none. In other words, such a 
system will behave very poorly with respect to load-balancing, and will give no good 
hint on the parallelization efficiency under normal (homogeneous) conditions. 

We hence abandoned this state point and instead systematically studied the weaker 
coupling Ib = 5, at first restricting the particle number to A^ = 500. Furthermore, we 
slightly increased the particle friction coefficient to C = 1-5- AH other parameters 
(density, temperature) were left unchanged. This is well in the homogeneous phase, 
and the pair correlation functions show much less structure, see Figs. E] and El It also 
turns out that here a larger MEMD lattice spacing a = 0.88 is sufficient to reasonably 
approximate the P^M result. 
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Figure 3. Pair correlation function of like charges at density p = 0.07 and Bjerrum 
length Ib = 5, comparing data obtained with P'^M with those from MEMD for different 
lattice spacings. 
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Figure 4. Pair correlation function of unlike charges at density p ■ 
length Ib ~ 5, again comparing P'^M with MEMD. 
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Figure 5. Internal electrostatic energy for P'^M, as a fmiction of the accuracy 
parameter AF/F. 



Our next aim is to compare the efficiency of P'^M and MEMD. In order to do this in a 
meaningful way, it is necessary to make sure that (i) both methods use parameters which 
yield roughly the same accuracy in the representation of the electrostatic interaction, 
and that (ii) both methods use parameters for which the results are obtained most 
quickly, within the accuracy constraint. The (thermally averaged) electrostatic energy 
[/ is a variable which, on the one hand, is easy to evaluate, and, on the other hand, 
reasonably sensitive to the long-range correlations between the particles. We therefore 
used this observable for calibrating the accuracy of the simulations. For P^M, we 
therefore calculated f/ as a function of the accuracy parameter AF/F. The results 
are shown in Fig. The error bars were obtained as statistical error bars, using the 
block average method [221 • Fi'om these results, one sees that an accuracy parameter of 
AF/F = 3.7 X 10^^ is good enough. This corresponds to the following P^M parameters: 
Mesh size 16^, third-order charge assignment, real-space cutoff 4.4, a = 0.43. 

At this point, it is necessary to comment on the evaluation of U in MEMD. The 
electric and magnetic field energies are given by (eo/2) / d^rE"^ and (l/(2eoc^)) / d^fH"^, 
respectively. Both types of fields have one longitudinal and two transversal degrees of 
freedom per lattice site. The longitudinal magnetic degrees of freedom are however not 
excited, due toV -H = 0. Furthermore, the magnetic part of the Hamiltonian is strictly 
quadratic, and the equipartition theorem can be applied. For this reason, the thermally 
averaged magnetic field energy is just given by MhsT, where M is the number of lattice 
sites. We have checked this relation, and found good agreement, except for a deviation 
of a few percent, which decreases with the time step, and must hence be attributed to 
discretization errors — the exact Boltzmann distribution is only generated in the limit 
of vanishing time step. This finding is a strong support of our belief that, except for 
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Figure 6. Internal electrostatic energy for MEMD, as a function of the lattice spacing 
a. 



the obvious conservation laws for the longitudinal fields, there are no further hidden 
conserved quantities in the system, and a thermostatting of the magnetic field is not 
necessary. We can also apply the equipartition theorem to the transversal part of the 
electric field energy, and hence the thermally averaged Coulomb energy is given by 

f/ = I / d'r (e^) - MkBT - (Useif) , (62) 

where Useif is the self-energy, as discussed in Sec. El We have measured f/ as a 
function of the lattice spacing a, using this recipe. While the results were in reasonable 
agreement with the P^M results, we did not find convergence for a — >■ 0. Rather, it 
seems that U diverges for small a (apparently like 1/a, though the data are not precise 
enough to be sure). It turns out that this divergence is reduced by reducing the time 
step, i. e. it is again an effect of discretization errors. Our explanation is that the 
cancellation of the self-energy is not perfect, because the subtraction term assumes 
the exact Coulomb lattice Green's function, while the simulation produces an effective 
lattice Green's function, which is slightly distorted by discretization errors. This effect 
is crucial for the energy calculation (and probably also for the evaluation of the pressure, 
and related quantities), but not for the particle configurations, which are stabilized by 
the repulsive LJ interactions. Probably that problem must be solved by combining 
MEMD with a Monte Carlo procedure, which can enforce strict detailed balance, and 
thus produce the exact Boltzmann distribution. For the moment, we just solved the 
problem by taking the particle configurations produced by MEMD, and using accurate 
P^M for evaluating U. The results, which now do converge, and give good agreement 
with the P^M data, are shown in Fig. IHl Again, we attribute the small remaining 
systematic deviation to discretization errors. Taking all these considerations together. 
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Figure 7. Particle difFusion constant as a function of the speed of light c. 



we took the lattice constant a = 0.88 as a value which produces sufficiently accurate 
data, consistent with our findings from the pair correlation function. 

Having fixed the parameters for accuracy, we now turn to optimizing the speed. In 
P^M, this was already done by the automatic routine (see above). In MEMD, we still 
have the speed of light c at our disposal. This parameter only infiuences the dynamics 
of the system, but not the statics. We actually checked that U does not systematicUy 
vary with c within our error bars — such a dependence could still be possible as a 
result of discretization errors. Furthermore, the CPU time necessary for one update 
step does not depend on c. Hence one would like to take a value of c for which the 
configurations decorrelate particularly quickly. In principle, each observable and its 
time autocorrelation function would have to be considered separately [221 • This is of 
course impractical, and hence we have taken the simpler criterion that the diffusion 
constant D of the particles, obtained from their mean square displacement, should 
be maximized. The corresponding data are shown in Fig. [3 One sees that D first 
increases as a function of c, but then saturates at a value which is in good agreement 
with the P^M value, except for some discretization errors. This nicely confirms the 
expectation that the dynamic properties should converge to electrostatic behavior for 
c — i> oo. We thus have the very favorable situation that computational efficiency for 
the statics and reasonable reproduction of the dynamics are not mutually exclusive, but 
rather congruent. However, this does not mean that one should just simply take a huge 
c value. Rather, c has to be small enough such that the Courant stability criterion 
|2n|, c <^ a/h, is still satisfied. For our parameters, a/h = 88, and hence c should be 
significantly smaller (in fact, our program crashed for c = 55). Therefore, we use the 
value c = 20. 



Coulomb Interactions via local MD 
1.00 



0.90 



21 



o 
CO 



0.80 



0.70 



0.60 



0.50 




10 15 20 25 

Number of processors 



Figure 8. Scalability factors s as a function of the number of processors, for both 
NEMD and P^M. For further details, see text. 



For the calibrated and optimized parameters, it makes sense to look at the speed 
in terms of CPU time. We ran the system on an AMD Athlon MP 2000+ processor 
for 2000 MD time steps. Since the diffusion constants for both methods are essentially 
identical, we do not need to take into account different rates of decorrelation. For P^M, 
the run used 17 seconds of CPU time, while for MEMD 16 seconds were needed. This 
shows that MEMD for a system of this density is a comptetitive alternative to P^M. 

Furthermore, we studied the scalability of our parallel programs at this state 
point. To this end, we systematically increased the particle number and the number 
of processors such that each processor keeps A^ = 4000 particles on average. The 
simulations were run on an IBM Regatta H server, where 10^ time steps were used for 
equilibration, and another 10"^ steps for measuring the CPU time. The P^M parameters 
(lattice constant of the mesh, real-space cutoff, a) were left unchanged throughout, 
since (i) one should expect that these values are reasonably close to the optimum also 
for larger systems, according to Refs. fl^ |^, and (ii) the automated optimization 
routine does not work very well for a very large number of particles. For a single 
processor, the timings were 759 seconds for P^M and 398 seconds for MEMD. We do 
not know for sure why the relative efficiency is so much different compared to the AMD 
processor, but speculate this might have to do with more efficient handling of memory 
for the Regatta architecture (note that MEMD is quite memory-intensive, due to all 
the various field variables). Figure |H1 presents the scalabity factors as a function of the 
number of processors, for both methods. The scalability factor s is defined as 



r{NpNy 



(63) 
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where N is the number of particles on a single processor, A''^ the number of processors, 
and r the total CPU time for a given number of steps. Note that for P^M we can produce 
reasonable data only for a processor number which is a power of eight, due to the Fast 
Fourier Transform in each spatial direction. For eight processors, we find that the 
scalability of MEMD is slightly better than for P^M. Furthermore, for MEMD we find 
reasonably acceptable (though not excellent) scalabihty behavior up to 32 processors. 

6. Conclusions 

MEMD is rather easy to implement and to parallelize. The numerical results, though 
being far from conclusive yet, seem to indicate that the algorithm is a competitive 
alternative to existing schemes for sufficiently dense systems. However, in electrostatic 
problems one often goes to much smaller densities. If we would apply the present MEMD 
method to such a dilute system, the number of grid points would become overwhelmingly 
large. P^M does not have this problem; due to the split-up of the work between real 
space and Fourier space it is possible to keep the number of grid points reasonably 
small. It is therefore clear that MEMD for such systems can only be competitive if it 
is also possible to use a reasonably coarse grid. We believe that this might be possible 
by introducing Yukawa subtraction combined with our Green's function subtraction for 
both the unscreened and the screened interaction. This further optimization of the 
method is left for future investigation. Another problem which needs to be addressed 
is the consistent handling of the discretization errors in calculating the energy, and 
related quantities — as we have seen, these interact in a very unfavorable way with 
the self-energy problems. We believe that this can be solved by combining MEMD 
with Monte Carlo, such that the Boltzmann distribution is reproduced exactly, and 
the potential of mean force is known exactly. Moreover, the dynamic properties of the 
algorithm have to be studied in more detail. In particular, it is necessary to investigate 
the accuracy of momentum conservation, and how this depends on the lattice spacing 
and the speed of light. This latter question is particularly important when considering 
applications which aim at dynamic properties, like, e. g. the dynamic behavior of 
charged colloidal suspensions. Much remains to be done, but the existing results are 
reasonably encouraging. 
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Appendix A. Details about Discretization 

A particularly useful spatial discretization scheme |S1 EI works as follows: The charges 
are interpolated onto the vertices r of a simple-cubic lattice with lattice spacing a. If 
the charge e^ is located at position fi in continuous space, then some nearby sites f 
are assigned some partial charges qi{f) = eis{f,ri) (s denoting a "smearing" function) 
such that J2r1i{^ = ^i or J2rS{r,ri) = 1. The total charge on site r is the sum of the 
contributions from all particles, g(r*) = J^iQii'^^ and the charge density is written as 
p{r) = a~^q{r). Different choices for s are possible; we have chosen linear interpolation 
to the eight vertices which form the cube in which the particle resides: 



s{r,ri) = (l \x -Xi\] 



X {l-^\y-y^\) (A.l) 

X (1 \z — Zi 



here x, y and z denote the lattice coordinates of the vertices. 

Now, the vector fields j, A and E are put on the links which connect the vertices, in 
such a way that they are aligned with the links. For instance, a link oriented along the 
X axis would contain a variable E which is the electric field at the position of the link, 
and which is positive if E points into the +x direction, while it is negative if it points in 
the —X direction. The divergence of such fields is put onto the sites, such that one just 
takes the differences of the field values associated with those six links which are directly 
connected to the site. Conversely, the curl of such fields is put onto the plaquettes by 
taking differences from the four field values which encircle the plaquette. The result is 
a vector perpendicular to the plaquette; positive (negative) field values are associated 
with a vector pointing in the +x {—x) direction (for the case that the plaquette is 
perpendicular to the x axis). Obviously, the fields and H must be such plaquette 
variables. The curl of plaquette fields is put onto the links, by taking differences from 
the four plaquettes adjacent to that link. Finally, the divergence of a plaquette field 
is put into the center of the cubes, by taking differences from the six plaquettes which 
enclose the cube. With these definitions it is easy to see that the divergence of a 
curl vanishes identically, as in the continuum, both for link and for plaquette fields. 
Furthermore, we can define the gradient of a scalar field, the latter being on the sites. 
The result is put on the links and obtained by just taking the difference between the 
field values on the adjacent sites. With this definition, one sees that the curl of a 
gradient vanishes, too. These identities are extremely important, since they allow us to 
decompose fields uniquely into longitudinal and transversal components, and to apply 
standard procedures of vectorial calculus also on the lattice. 

The particle motion generates currents on the surrounding links. We again use a 
linear interpolation scheme for j, where the current is distributed onto the twelve links 



Coulomb Interactions via local MD 24 

which surround the cube in which the particle is: For a hnk / oriented in the x direction 
the current contribution from particle i is 

j{l) = a-'^avi^ {wi{l) + W2{1)) (A.2) 

where Vi^ is the x component of the particle's velocity, while wi and W2 are the charge 
weight factors of the two sites which are connected by the link. For the y and z direction, 
the analogous procedure is applied. It is easy to see that the space-discretized continuity 
equation holds exactly. Similarly, the discretized version of the fourth Maxwell equation 
(Eq. E]) implies that the discretized version of Gauss' law (Eq. Q) holds exactly for all 
times if it holds at time t = 0. The discretization scheme is therefore suitable to keep 
the system on the constraint surface. 

Apart from interpolating the charges and currents onto the lattice, we also need to 
interpolate the electric field onto the particles in order to calculate the electric force. 
Here we use the same scheme as for the current (Eq. IA.2J1 . i. e. the field in, e. g., 
the X direction is obtained by summing the fields from the four surrounding links in x 
direction, weighted by the sum of the two charge weight factors of the sites connected 
by that link. 

Using this scheme, the whole theory of Sec. Elis consistently discretized. The system 
is initialized by putting particles into the simulation cell (which has periodic boundary 
conditions), assigning velocities and charges to them (of course, the overall system is 
neutral), and calculating the electrostatic electric field as follows: First, we exploit the 
fact that the solution of Eq. ^ is trivial in one dimension. This allows us to find a 
simple solution by just treating the spatial dimensions recursively: First, the field in z 
direction is calculated by taking into account the differences between the mean charges 
of planes perpendicular to z. Within the planes, we then take into account charge 
differences between lines (after subtracting the average plane charge) to obtain the field 
in y direction. Finally the field in x direction is obtained from the charge differences on 
the sites within a line. This solution of course violates V x E = 0. In order to bring 
the system onto the BOS, we iteratively relax the 9 field on the plaquettes until the 
electrostatic field energy is minimal (cf. Eqs. 1^ fT^. For a single plaquette, this can 
be done in a single step. For the overall system, we use a checkerboard decomposition 
which allows easy parallelization. The field A is initialized as zero. Then Eqs. IHT^HUI 
are integrated in time. 

Appendix B. Integrator 

Ideally, one would like to run MEMD via an integrator which leaves the phase-space 
volume invariant and is time-reversible, such as the Verlet algorithm in standard MD 
pij . Since the equations of motion (even in the lattice-discretized case) have these 
properties, it is indeed possible to construct such a scheme. Disregarding Yukawa 
subtraction for the time being, an analog to the Verlet algorithm for MEMD would 
be the following integrator for Eqs. IHTHIUl based upon a time step h: 
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(iv 

(v 
(vi 

fvii 



Update the particle momenta by half a time step. 

Update the A field by half a time step. 

Update the particle positions by half a time step. 

Update the electric field by a full time step. 

Update the particle positions by half a time step. 

Update the A field by half a time step. 

Update the particle momenta by half a time step. 



Here, "update" means the simple Euler rule x{t + h) = x{t) + x{t)h. The time 
consuming part (update of the particle momenta, update of the electric field) is arranged 
in such a way that only one "force calculation" per time step is necessary. This scheme 
does conserve the phase-space volume and is time-reversible, however, it suffers a 
severe disadvantage: The update of the electric field (step 4) is based upon a particle 
configuration (in real space and velocity space) which has so far only progressed by half 
a time step. As a consequence. Gauss' law is not satisfied within machine accuracy, but 
rather only within the accuracy of the time discretization (to satisfy it exactly would 
require to know the current at the end of the time step, too). This is very undesirable, 
and hence we have adopted the elegant solution which was found by Rottler and Maggs 
PH] and allows to conserve both time-reversibility and phase-space volume conservation, 
while keeping the system strictly on the CS: 
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(vii 

(viii 

(ix 
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(xii 

(xiii 

(xiv 
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Update the particle momenta by half a time step. 
Update the A field by half a time step. 

Update the particle positions in x direction by half a time step. 
Update the electric field in x direction by half a time step. 
Update the particle positions in y direction by half a time step. 
Update the electric field in y direction by half a time step. 
Update the particle positions in z direction by half a time step. 
Update the electric field in z direction by a full time step. 
Update the particle positions in z direction by half a time step. 
Update the electric field in y direction by half a time step. 
Update the particle positions in y direction by half a time step. 
Update the electric field in x direction by half a time step. 
Update the particle positions in x direction by half a time step. 
Update the A field by half a time step. 
Update the particle momenta by half a time step. 

We have added a Langevin thermostat to the particles: 

|p. = -|| + e.^(r-;)-^pl + /; (B.l) 
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where ( is the particle friction constant, and fi is a random force satisfying the standard 
fluctuation-dissipation theorem: 

fmfjit')) = 2CkBT5,,5ap6it - t'), (B.2) 

where a and f3 denote Cartesian indices. This puts the system into the canonical 
ensemble. For large systems, one can rely on the equivalence of ensembles, and there is 
no fundamental statistical-mechanical need for such a thermostat — it is just a matter of 
technical convenience: Usually a Langevin thermostat tends to stabilize the simulation 
due to its inherent feedback mechanism, such that larger time steps are feasible. Rottler 
and Maggs (TU] also add a Langevin thermostat to the magnetic field; we have not done 
this. It should be noted that such thermostatted dynamics violates time reversibility 
and phase-space volume conservation anyways. 
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